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I. INTRODUCTION 



For a structural analysis of any non-trivial elasto-static 
problem by the finite element method, a large amount of data must 
be provided to the computer system which performs the calculations. 

The required information includes the arrangement of nodes within 
each element, the coordinates of all nodes, and the applied loads. 
After the numerical values of this data have been determined, which 
is a large problem in itself, the information must be supplied to 
the computer system. This is normally accomplished by punching the 
appropriate data on computer cards. The typing of these cards is 
also a difficult and time-consuming job. More important, however, 
is the ever-present possibility of human errors during the performance 
of these tasks. If an error is not detected prior to submitting 
the data for processing, expensive computer time would be wasted. 

With the above considerations in mind, it is readily apparent 
that a computer system, which would generate most of the data required 
by the primary computer system for a complete structural analysis, 
would be very useful. In addition, this mesh generator should 
require as little input data as possible. 

References 5, 8, and 10 contain reports of efforts which have 
been directed at coding computer systems which will aid the analyst 
in the interpretation of output data. A few computer systems, which 
are reported in Refs. 1, 2, and 9, have been coded to automatically 
generate input data. However, these input data generation systems 
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are either for types of finite elements other than isoparametric 



elements or for tvo-dimens iona 1 eLements, 

Two input data generator systems for three-dimensional i 
parametric finite elements will be presented in this thesis, 
system is for quadratic elements and the. second system is for 
elements . 



o- 

One 

cubic 
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II. ISOPARAMETRIC FINITE ELEMENTS 



The isoparametric concept will be discussed in this chapter in 
sufficient detail to establish only the procedures for constructing 
isoparametric elements and the validity of relationships used to 
obtain the consistent gravity and pressure loads. The readers are 
referred co Refs. 4 and 7 for further details of the isoparametric 
concept . 



A. DISPLACEMENT FUNCTIONS 



The displacements of nodes in any one element are obtained by 
using shape functions which are defined in a system of dimensionless 
coordinates (§,T|.Q that are unique for that element. The coordinates 
range from minus one to plus one. 
u(x, y, z) = £ N i (5,Tl,0 

v(x, y, z) = ?N. (5,11,0 v. 

111 ( 1 ) 



w(x, y, z) =SN.(§,11,0 w. 

i 1 1 



where 1L (£,!],£) are the shape functions and u_^, v^, w^ are the dis- 
placement components of node "i" in a global reference system (x , y, 

where ,f i" assumes values from 1 to the total number of nodes in the 

► 

element. This global reference system is then related to the dimen- 
sionless coordinates by the same shape functions as follows: 
x(5,71,0 = s N. (5,11,0 X t 

y(5,H,0 « Sn (5,11,0 y t 

1 ( 2 ) 



z(5,H,C> ■= 5 n. (|,il,0 

i i 1 



z) , 
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where x^ , y , z^ are the coordinates of node "i ,f in the global 
reference system [4]. 

The shape functions referred to above have been obtained from 
[7]. They are listed in section 2.C. 

B. TRANSFORMATIONS 



In order to obtain the consistent gravity and pressure loads, 
the volumetric increment, dxdydz, must be found in terms of the 
dimensionless coordinates. That operation requires the Jacobian 
matrix [4], 



[j] = 



3x/S5 dy/3§ Sz/S§ 

ax/si] ay/aii az/aii 

ax/ac ay/ac az/ac 



(3) 



The element of volume is then transformed as follows: 

dV = det[j]d§ dTj dC ' (4) 

C. SHAPE FUNCTIONS 

In this listing of the shape functions, attention is directed 
to the tri-quadratic and tri-cubic elements, for which the mesh 
generator programs were written. 

1 . Quadratic Elements 

See Figure 1 for identification of the nodes. 

Corner nodes: 1, 3, 5, 7, 13, 15, 17, and 19. 

n. = (i/s) (i + c o ) (i + *n o ) (i + e 0 ) <s 0 + ti o + c o - 2 ) (5) 

where = 5^5 and ^ is j* 1; similarly for Tj^ and 
Mid-side nodes: 2, 6, 14, and 18. 

N. = (1/4) (l - S 2 )(l + n o )(l + C D ) (6) 
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( 7 ) 



Mid-side nodes: 4, 8, 16, and 20, 

N. = (1/4) (l - Tj 2 ) (1 + § o )(l + C c ) 

Mid-side nodes: 9, 10, 11, and 12. 

N. = (1/4) (l - C 2 )(l + S 0 )U + \) (8) 

2 . Cubic Elements 

See Figure 2 for identification of nodes. 

Corner nodes : 1, 4, 7, 10, 21, 24, 27, and 30. 

N. = (1/64) (l + § o )(l + T| q ) (1 + C 0 )[9(§ 2 + Tl 2 + C 2 ) - 19] (9) 
Mid-side nodes: 2, 3, 8, 9, 22, 23, 28, and 29. 

N. = (9/64) (1 - l 2 ) (1 + 9§ o ) (l + Tl o ) (1 + C o ) (10) 

where = t 1/3, T) = * l, and C, = * 1. 

Mid-side nodes: 5 , 6, 11, 12, 25, 26, 31, and 32. 

N. = (9/64) (1 - T] 2 )(l + 9T] )(1 + § )(1 + C ) (11) 

1 GOO 

where Tj. = .1" 1/3, = - 1, Q. = + 1. 

i l l 

Mid-side nodes: 13, 14, 15, 16, 17, 18, 19, and 20. 

N. = (9/64) (1 - C 2 )(l + % )(1 + § )(1 + 'll ) (12) 

1 o o o 

where C- = - 1/3, I = t 1, and 11. = - 1. 

Note that the convention for numbering of nodes within any 
element is to begin at the location where § = T] = £ = 1 and to 
consecutively number the nodes, proceeding in a counter-clockwise 
direction around the £-axis (see Figures 1 and 2). 
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15 14 13 




Figure 1. Quadratic Element 
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24 23 22 21 




Figure 2. Cubic Element 



III. DISCUSSION OF COMPUTER PROGRAMS 



The computer programs presented in this thesis were written to 
support a computer system called TRISOP. TRISOP was coded by Professor 
G. Cantin of the U. S. Naval Postgraduate School, Monterey, California. 
TRISOP performs a structural analysis of three-dimensional elasto- 
static problems using isoparametric finite elements. Information 
required by TRISOP for any problem includes general mesh parameters, 
title cards, format statements, and material properties; Young’s 
modulus of elasticity, Poisson’s ratio, and the coefficient of thermal 
expansion. Fifteen cards are needed to supply this data to TRISOP. 

The bulk of the data required by TRISOP is the element connectivity, 
structure geometry, loading, and boundary conditions. The boundary 
conditions will vary with each- problem; therefore, no effort has 
been made to automate this process. 

The mesh generating computer programs discussed in this chapter 
compute the remainder of the data required by TRISOP: element 

connectivity, coordinates of nodal points, consistent gravity load 
vector, and consistent pressure load matrix. Additional information 
which is helpful in preparing an input deck of cards for TRISOP is . 
also computed. One of these programs generates data for quadratic 
elements and the second program generates data for cubic elements. 

These programs will be called QUAMEG and CUMEG, 

A. ELEMENT CONNECTIVITY 

A structure of any shape can be visualized as a cubical piece 
of plastic material which can be stretched, compressed, and otherwise 
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appropriately deformed in such a way that it will assume the geometry 
of the structure. The element connectivity is determined when the 
structure is in this initial abstract cubical shape. The dimensionless 
coordinates of this abstract structure are §*, T] / , (J , 

Numbering of nodes is begun in a face which will have the 
fewest nodes. In that face, numbering proceeds sequentially along 
the edge having the smallest number of nodes as shown in Figure 3. 

This convention reduces the size of the half-bandwidth of the system, 
which directly affects the space and time required by the computer 
to solve the problem. The coordinate system is positioned such 
that node one is at %> ' = 0, T] 7 = 1 , and Q* = 1. An example is 
discussed in Appendix A. 

For this mesh, the connectivities of the elements follow simple 
but tedious rules that were coded directly. 

B. SOLUTION TIME AND SPACE REQUIREMENTS 

The time required to solve the system of equations in the problem 
and the space required by the computer to store all data are important 
information. They must be determined before the problem is submitted 
for processing. The calculation of this information requires no 
additional input data and it is therefore included in QUAMEG and CUMEG 
to reduce the effort required of the programmer. 

The time required to solve the equations is computed using 
an empirical formula which was constructed by Cantin [4]. This 
calculation has been included in the coding of QUAMEG and CUMEG. 

The space required to store the various bits of information 
for a problem depends upon the total number of records and the 
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Figure 3 . Abstract Structure 



length of each record. These will vary with different types of 
computers. Simple formulas which generate the correct results for 
the IBM 360/67 computer used at the U. S, Naval Postgraduate School 
are employed in the coding of QUAMEG and CUMEG. 



C. NODES ON BOUNDARY SEGMENTS 



Before the coordinates of nodes which are within the boundaries 
of a structure can be computed, the coordinates of ail corner nodes 
on the edges of the bounding surfaces must be specified. This problem 
is simplified by the method discussed below. 

The term "boundary segment 11 must first be defined for the purpose 
of this discussion. A segment must satisfy a few simple conditions 
which depend upon the orientation of the proposed segment in the 



structure and upon 



he type of coordinate system used to describe 



the geometry of the structure,. QUAMEG and CUMEG are coded to allow 
the use of either rectangular or cylindrical coordinate systems. A 
combination of the two coordinate systems in any one problem is not 
permitted . 

1 . Rectangular Coordinates 

When rectangular coordinates are used, a segment is any 
straight line between two corner nodes of the mesh. The element 
divisions along this line must be uniform. There are no other restric- 
tions imposed upon a segment when rectangular coordinates are used. 

2 . Cylindrical Coordinates 

When cylindrical coordinates are used, the radius and reference 
angle are taken in planes parallel to the x-y plane. In this case 
a segment can be an arc of circle centered on z and in a plane 
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parallel to plane x-y, or a straight line in a plane containing 
the z-axis. Again, element divisions along these lines must be 
uni form , 

In order to determine the coordinates of nodes along a boundary 
segment, the coordinates of the corner nodes at the ends of the 
segment must be provided in the input data. In addition, the location 
of the segment within the mesh must be specified. The method of 
specifying the location of a segment is described in sections 4,R,3 
and 4,B,4, 

In Figure 4, the vector AB is easily obtained from the coordinates 
of nodes A and B. After calculating the values of increments of 
distances along x, y, and z, the intermediate nodes are easily 
obtained . 

If the bounding surfaces of a structure are of an irregular 
configuration such that portions of the edges of that surface cannot 
be specified by segments, the coordinates of all nodes along this 
irregular portion of the bounding surface must be provided in the 
input da ta , 

D, CORNER NODES 

Once the coordinates of the nodes on all of the boundary segments 
are established, the coordinates of corner nodes on the bounding 
surfaces and on interior slices are calculated using a method described 
below. The bounding surfaces are considered in the following order 
(Figure 3): face 1 ((/-l) , face 2 (C /= 0), face 3 (§ /= 0), face 4 (§ /s =l), 
face 5 (T| ; = 1 ) , face 6 (T| 7 =0) , Slices which are each of constant Q* 
are then considered until the coordinates of all corner nodes are 
calculated , 
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Figure 4. Boundary Segment 
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The method used for calculating the coordinates of corner nodes 
within the edges of a bounding surface or slice can best be understood 
with the example illustrated in Figure 5. In this example, the 
surface ABCD is the real surface for which a mesh is desired. The 
mesh is first constructed in the plane as shown. 

At this time, the coordinates of all the corner nodes on the 
boundary are known. If the coordinates of the interior nodes are 
required to satisfy the Laplacian equation: 



V (x) = 0 
V 2 (y) = 0 
V 2 (z) = 0 

2 2 2 2 2 
where V = 9 /Scv + 9 /9j3 



(13) 



a set of acceptable interior nodes will be obtained. The dummy 
variables (o', 8) are replaced by appropriate combinations of § / , Ij* , 
and C ? / , depending upon the location of the face being considered. 
For this computation, the standard Laplacian molecule of finite 
difference theory is used and the final coordinates are obtained 
by relaxation. Since the coordinates of the nodes within the edges 
of a face are arbitrary, the relaxation process can be stopped at 
any time, depending upon the degree of approximation desired. In 
this code, sixty iterations are used. 



E. MID-SIDE NODES 

After the coordinates of all corner nodes have been established, 
the coordinates of the mid-side nodes are computed by averaging 
the coordinates of adjacent corner nodes. Although the mid-side 
nodes need not be located exactly mid-way between the corner nodes, 
it is the simplest means of determining their coordinates. 
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Figure 5 . 



Bounding Surface 
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F. CONSISTENT GRAVITY LOAD VECTOR 

The z-axis has been chosen as the line o£ action of the gravity 
load. The work done by the force due to gravitational acceleration 

on a. body is 

n (14) 

U = \ v pgwdV 

whe re V is the work done; p is the mass density of the body; g is 
gravitational acceleration; w is the displacement of the body in 
the z-direction; and dV is the element of volume. 

The work done by concentrated forces acting at each node in 

the structure mesh is, in matrix notation. 

[ 1 (15) 

U - <V>{w.j> 

where (V) is the consistent gravity load vector and ^ are e 

respective node displacements. 

From Chapter 2, recall the following two expressions: 

w(x,y,z) = ? N i (5,Tl,C)^ w i | 
dV = det[j] d^dTjdC. 

When these expressions are substituted into equation 14 and the 
result identified to equation 15, the consistent gravity load vecto 

is obtained, 

<V> - PS \ l ( t iet[J3<N )d5dHdC (16) 

-l-l-l 

In this form, equation 16 is ideally suited for solution by 
Guassian cubature. In QUAMEG, four Gauss points are used in each 
of the three directions, for a total of 64 points within each element. 
Five Gauss points are used in CUMEG. The weight factors and coor- 
dinate values are obtained from Ref. 11. The Jacobian matrix is 
determined by using three function statements and the shape function 
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are determined with two function statements in Subroutine GRAP 
(Appendix C) . These function statements are permuted to obtain the 
partial derivatives of the shape functions and the values of the 
shape functions for all nodes in the mesh. 



G. CONSISTENT PRESSURE LOAD MATRIX 



The consistent pressure load matrix for a uniformly distributed 
pressure acting on any face of an element is determined in a manner 
similar to that used for the gravity load. 

The work done by the uniformly distributed pressure is: 
u = ^ A (-pn dA) - GI) (17) 



where n is a unit normal vector (positive outward) and H is the 
displacement vector at the surface (see Figure 6). The displacement 
vector u can also be written as: 

u = u(x,y,z)T + v(x,y,z) j + w(x,y,z)k (18) 

The vector ndA is found from differential geometry: 

ndA = (37/dcOX 07/36) dadg (19) 



where r = xi + yj 4- zk is a position vector of a point on the pres- 
surized surface and a and f3 are dummy variables, standing for either 

I, “H, or Q. 

The terms 37/3(2 and 37/38 are simply two rows -of the Jacobian 
matrix evaluated on the surface. After making appropriate substi- 
tutions and performing a few matrix manipulations, the expression 
reduces to: 



U 




<f(a,P)> 




dcfdg 



( 20 ) 



where 



• 2 * 



is the vector of nodal coordinates for the eight nodes 



on the face and f^,^) is a row vector of functions of a and 0. 
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Figure 6, Pressurized Element 



The integration is performed by Gaussian quadrature and the 
resultant vector is exhibited in convenient matrix form. Five Gauss 
points are used in both QUAMEG and CUMEG. 

The Jacobian matrix and shape functions are determined in 
Subroutine GRAP by using the same function statements as are used 
in the calculation of the consistent gravity load vector. The face 
number of the element on which the pressure is acting must be included 
in the input data. With this information available, a "computed 
go to" statement is used in the code to determine which variables 
(5, T], O will replace ct and (3. Referring to -Figure 6, the faces 
are located as described below. 



Face 


Location 


21 


3 


Face 


Location 


a 


3 


i 


§ = +i 


T] 


C 


4 


§ = -i 


c 


71 


2 


T| = +1 


c 


§ 


5 


T1 “ -l 


i 


£ 


3 


C - +1 


§ 


n 


6 


C = -1 


71 


§ 



H. PLOT OF STRUCTURE MESH 

Once the coordinates of all nodes in the mesh have been established, 
the mesh must be inspected to determine if it is acceptable. Without 
a plot of the coordinates, this task is nearly impossible. Plotting 
the mesh by hand is prohibitively time-consuming. Subroutine GRID 
has been included in QUAMEG and CUMEG to generate an off-line printer 
plot of the mesh. The parameters required for IBM 360 source library 
program DRAW are determined in GRID. DRAW is called to perform the 
actual plotting of the mesh. GRID generates a two-dimensional 
plot of the mesh on slices at constant values of Q 1 . For the example 
shown in Figure 3, faces 1 and 2 and five slices at constant values 
of Q / would be plotted. 
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The user ran specify that the scaling be performed automatically 
by DRAW, However, that choice could result in the x-scale being not 
equal to the y-scale. The plot of the mesh would be distorted in 
that case. The scaling is done in GRID to insure that the scales 
will be equal and that the smallest possible scale is used. 

The plotting space available with DRAW is fifteen inches in 
the y-direction and nine inches in the x-direction. Recall that 
the convention for numbering the nodes usually results in having 
the narrowest side of the structure along the y-axis (section 3. A), 
Therefore, GRID interchanges the x- and y-coordinates of the structure, 
for plotting purposes only, to take advantage of the larger plotting 
space in the y-direction. This further reduces the scale of the 

i 

plot. 
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IV. INPUT DATA PREPARATION 



This chapter is written as a self-contained unit giving specific 
information necessary for the preparation of input data cards. The 
chapter may be removed and used as a guideline for input data preparation 
after a thorough knowledge of the computer programs has been obtained. 



A. VARIABLES 



The variables and formulas defined below are included in the 
code and should be used as an aid for determining the arrangement 
of nodes within the mesh when preparing input data. 

Variable Meaning 

L. NEL Total number of elements. 



NEL = (NX) (NY) (NZ) , where NX, NY and NZ are the number of elements 
along the £ ,T| , and £ directions (21) 

2. NPS Total number of nodes on any face at constant 

£ , where a face is a slice which intersects 
both corner nodes and mid-side nodes. 



NPS « 3 (NX) (NY) + 2 (NX + NY) + 1 (QUAMEG) (22) 

NPS - 5 (NX) (NY) + 3 (NX + NY) +1 (CUMEG) - (23) 



3 . NPM 



Total number of mid-side nodes between any two 
faces. 



NPM = (NX + 1) (NY + 1) (QUAMEG) 

NPM = 2 (NX + 1) (NY + 1) (CUMEG) 



(24) 

(25) 



4 . NPL 

NPL = NPM + NPS 

5 . NUMNP 



Total number of nodes 
layer consists of one 
intersecting mid-side 
and an adjacent face. 



Total number of nodes 



in one layer, where a 
face plus the slice 
nodes between that face 

(26) 

in the mesh. 



NUMNP = (NZ) (NPL) + NPS 



(27) 
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B# DATA CARDS 

!. Element Field °at a (7110,215) 



Number of cards per i 
Columns Variable 


l - 10 


NX 


ll - 20 


NY 


21 - 30 


NZ 


31 - 40 


NDX 


41 - 50 


NDY 


51 - 60 


NDZ 


61 - 70 


NPDP 


71 - 75 


KCUD 



76 - 80 



NCARD 



Meaning 

Number of elements in the § -direction, 

Number of elements in the l] -direction, 

Number of elements in the C'-direction, 

Number of segments in the ^'-direction; 

Number of segments in the ^-direction; 

Number of segments in the C -direction, 

Totnl number of nodes^for ^ichUie coor 

tocation"o£ segment end points or other 
boundary nodes ; 

cal' coor dinates'are'used^ in the problem. 

coord ina t es ‘and^a' blank' or' ter t> L specifies 

rectangular coordinates; 

A positive integer indicates that punched 

hot desired. 



Alphameric information to identify the 
prob lem . 



2 . Problem T Hpnt if ication (10A8) 

Number of cards per problem - one 

Meaning 

Columns Variable 

l - 80 TITLE 

3 . Modal Point Coordd riates (UlO, 3F10.0) 

_ one for each NPDP point 

Number of cards per problem 

Meaning 

Columns Variable 

I Nodal point number; 

U - 20 COORD (I , l) x-coordinate (radius), 
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Columns 


Variable 


Meaning 


21 - 30 


COORD (1, 2) 


y-coordinate (angle) ; 


31 - AO 


COORD (I, 3) 


z-coordinate . 


4. Segment Divisions 


(5110) 


Number of 


cards per problem - one for each segment 


Columns 


Variable 


Meaning 


1-10 


i 


Lowest node number in the segment; 


11 - 20 


j 


Highest node number in the segment; 


21 - 30 


N1 


Segments in the §'- and £ '-directions : 
Number of elements above the segment 
(in the positive ^'-direction) ; Segments 
in the 7] '-direction: Number of elements 

within the segment; 


31 - A0 


Ml 


Segments in the T|'- and £ '-directions : 
Number of elements to the left of the 
segment (in the negative § -direction); 
Segments in the § -ditto lion ; Number 
of elements within the segment; 


Al - 50 


LI 


Segments in the - and 1) '-directions : 

Number of elements before the segment 
(in the positive £'-direction) ; Segment; 
in the £ '-direction : Number of elements 

within the segment. 


5. Load Data (2110, 


2F10.0) 


Number of 


cards per probl 


em - one 


Columns 


Variable 


Meaning 


1-10 


MAP 


A positive integer indicates that a 
plot of the structure mesh is desired. 
A blank or zero indicates that a plot 
is not desired; 


11 - 20 


NEFP 


Number of pressurized element faces; 


21 - 30 


SGZ 


Specific weight ; 


31 - A0 


UDP 


Uniformly distributed pressure. 
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6. Pressure Loading Data (2110) 

Number of cards per problem - one for each NEFP 



Columns 


Variable 


Meaning 


1-10 


NELP 


Element number on which a pressure exis 


11 - 20 


NFACE 


Face on which the pressure is applied. 
(See section 3.G and Figure 6) 


7. 


Plot Title Cards (6A8) 




Number of cards per problem - two per plot (NZ + 1 plots) 


Columns 


Variable 


Meaning 


1-48 


ITITLE 


Alphameric information to identify the 
plot; 


r- » 

1 

-O 

00 


ITITLE 


Alphameric information to identify the 
programmer and his computer center box 
number . 


8. 


Stop Trap 





The stop trap employed in the programs allows the execution 
of as many separate problems as may be desired with one submission 
of the program. The program will continue processing problems until 
it is informed that NX is less than zero. Therefore, the last card 
of a complete input data deck should have a negative integer in 

columns 1 - 10. 
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v. CONCLUSIONS 



The addition of QUAMEG and CUMEG to the tools available in the 
finite element method of analysis opens new doors for the analyst 
who uses isoparametric finite elements. Previously, the amount 
of input data required to analyze a structure of appreciable size 
was too overwhelming for the analyst to undertake the problem. 

QUAMEG and CUMEG make it practical for the analyst to solve such 
a problem. In addition, he can be certain, after an acceptable 
mesh is generated, that the output data generated by QUAMEG and 
CUMEG is free of error. 

For the example discussed in Appendix A, the number of input 
cards required by QUAMEG versus the number of output cards generated 
by QUAMEG for direct input to TRISOP is about six per cent. For 
the 1x8x8 mesh of the pinched cylinder problem discussed in 
Appendix B, this ratio is about 4 per cent. 

With a gravity load and pressure load acting on one cubic 
element, it was found that QUAMEG and CUMEG produced consistent 
load vectors which were accurate to at least eight decimal places. 

QUAMEG and CUMEG considerably reduce the effort required of 
the analyst to solve a problem by the finite element method, generate 
acceptable values for nodal point coordinates and compute accurate 
consistent gravity and pressure load vectors. 
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VI . RECOMMENDATIONS 



A. LOADING CONDITIONS 

The capabilities of QUAMEG and CUMEG should be increased by 
adding subroutines which will generate the consistent load vectors 
for a thermal gradient in the structure and for an initial displace- 
ment condition at one or more nodes. 

B. PLOT OF STRUCTURE MESH 

It would be preferred to obtain a more complete representation 
of the structure mesh than is now provided. Two courses of action 
are recommended: 

1. Expand Subroutine GRID in such a way that plots of faces 

in the and/or planes are also generated. 

2. Replace GRID with a subroutine which would generate a 
three-dimensional plot of the structure mesh. The feasibility 
of using the IBM 360 source library program PLOT PACKAGE, which 
is mentioned in Ref. 12, should be investigated if this method 
is attempted. 

This author feels that the results of recommendat ion 2 would 
be preferred to 1, although the coding may be more difficult. 
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APPENDIX A 

GENERATING A MESH WITH QUAMEG 



This appendix describes the preparation of input cards and the 
expected results for a sample problem using QUAMEG. This problem 
was chosen as an example to demonstrate the capabilities of QUAMEG 
and to more fully explain the procedures described in Chapter 4. 

The structure considered is a cylindrical solid with a star 
shaped void in the center. This void linearly converges to a point. 
Figure 7 shows a 45° sector of the structure. It is normally prefer- 
able to work with angles in the first quadrant. Therefore, the 
sector from zero . to forty- five degrees will be considered. 

The first step in the problem is to decide upon a mesh size. 

In this case, an 8 x 5 x 4 mesh was chosen. The amount of control 
to be exercised by the programmer to insure that the structure 
will be modeled accurately can be determined now. The coordinates 
of all nodes along the edges of the bounding surfaces of the structure 
must be available to the system before the remaining nodal coor- 
dinates can be determined. Since there are no irregularities on 
the boundaries between the front and back faces the coordinates 
of the nodes on the edges between these two faces can be computed 
by providing the necessary segment data. Therefore, only the nodes 
on the front and back faces need be considered. 

At this point, two rectangles of constant Q* are drawn with 
the proposed mesh superimposed upon them (see Figures 8 and 9) . 

These figures will serve as an aid for numbering the nodes on the 
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actual structure and for determining the locations of segments. 

After the nodes on face one are numbered, eauation 22 is solved to 
insure that an error has not been made in the numbering. The result 
of equation 22 is subtracted from that of equation 27 to determine 
the number of the first node on the back face. In this way, the 
programmer need not go through the time-consuming process of numbering 
all nodes in the mesh. The nodes on the back face are numbered and 
checked as before. Note that the 8 x 5 mesh is in the § -T| plane. 
Orienting the structure so that the 5x4 face of the mesh were in 
the § / -T| / plane would reduce the half-bandwidth of the system. 

However, cylindrical coordinates must be used in this problem and 
it would not be possible to model the structure accurately if it were 
oriented in that manner. Finally, a picture is drawn of the front 
and back faces of the structure with the desired mesh super imposed 
and the nodes numbered (Figures 10 and 11) . 

The identification of segments will now be determined. The 
segments will be considered in the same order as they are considered 
in the program. Referring to Figures 8 and 9, it is seen that the 
following lines are parallel to the V-axis: 1 “ U > 137 “ 147 ’ 

805 - 815, 941 - 951. Now referring to Figures 10 and 11, it is 
seen that all of these lines are of constant angle with a uniform 
distribution of elements within the lines. Therefore, these lines 
will be treated as segments and the coordinates of the nodes listed 

above will be included in the input data. 

By this method of specifying lines, the author does not intend 
to imply that the nodes along the line are numbered consecutively, 
starting with the lower number. For example, the nodes along line 
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60 - 111 in Figure 8 are: 60, 67, 77, 84, 94, 101, and 111, This 

method of specifying lines merely lists the numbers of nodes at the 
ends of a line. 

In order to control the departure angles of lines 79 - 77 , 

96 - 94, and 113 - 111, it is necessary to force line 52 - 62 nearer 
the x-axis than would result if an even distribution of nodes along 
line l - 137 were allowed. Therefore, an appropriate location for 
node 52 is chosen, the coordinates of nodes 52 and 62 are included 
in the input data, and line 52 - 62 is treated as a segment. Note 
that this line is of neither constant angle nor constant radius. 
However, since the line is not on the edge of a bounding surface, 
the fact that it will be curved will not affect the solution adversely. 
Similarly, the coordinates of nodes 856 and 866 will be included 
in the input data and that line will be a segment. 

The same rules govern the selection of segments in the -direc- 
tion. Referring to Figures 8 and 9, it is seen that the coordinates 
of nodes along lines 1 - 137 , 11 - 147 , 805 - 941 , and 815 - 951 
must be made available to the computer. In Figures 10 and 11, line 

I - 137 is a line of constant radius; however, the element divisions 
along this line are not uniform. Two segments must be specified: 
lines 1-52 and 52 - 137. The coordinates of the end nodes on these 
segments have been specified already and need not be repeated. Line 

II - 62 is of constant radius and is therefore one segment. Again, 
the coordinates of these nodes have been specified previously. Line 
62 - 113 is of constant angle, meeting the criterion for a segment. 

The coordinates of node 113 must be specified. Line 113 - 147 is 

of neither constant angle nor constant radius. Therefore, it cannot 
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t- and the coordinates of all nodes along this 
be treated as a segment and the coor 

u . r In order to further control 

line must be included in the input data. . 

the departure angle of the lines mentioned above, line 60 - 145 
Is also included as a segment. The coordinates of the end nodes 
o£ this segment have been computed in the problem already and do 

not need to be included in the input data. As with line 52 - 62. 

. £ ii ne 60 - 145 is not important. Similarly, 

the exact shape of Line 

„ns 856 856 - 941, and 864 - 949 are treated 

the back face, lines 805 856, 85b 

- A qk . 951 are at the 

as segments. Although the nodes along Un 

origin, care must be taken to insure that the .coordinates of these 
nodes are specified at appropriate angles in such a may that thi 
point mill have the same "shape” as the front face. The radial 
component of these coordinates is aero. This care is necessary 
hecause all computations in the problem are made «th cylindrical 
coordinates. Were the coordinates of the nodes along this line 

he distorted. The folding lines can also be treated as segments. 

Rfifi 917 The coordinates of nodes 923, 934, and 
815 - 866 and 866 - 9L/. me 

, • rtata for the same reason as vas 

940 must be included in the input data 

given for nodes 119, 130, and 136. 

To complete the process of specifying the coordinates of nodes 

along the structure edges, the following lines mill be treated as 

segments in the C'-directio„: 1-805,52-856 Ul - 941, 60 - 864, 

77 . 881, 94 - 898, 11-915, 11-815,62-866, 113 -912, and 
130 - 934. 

The procedure for determining the values of Nl , 
for the segment division data will not be discussed 
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felt that the explanation in Chapter 4, combined with the listing 
of the input cards for this problem, is sufficiently clear. 

The card containing the element field data is prepared after 
the above steps have been accomplished. Loading information and the 
request for punched cards should not be included in the input data 
until it has been determined that the mesh will be satisfactory. 

Approximately 45 seconds of computer time were required to solve 
this problem with QUAMEG . That time includes the 31 seconds required 
to compile the program. Figures 12, 13, 1.4, 15, and 16 are the plots 
of the structure mesh that was generated by QUAMEG. Table I is the 
listing of the required input cards described above. Figure 17 is 
the result of a 5 x 5 x 3 mesh generated with CUMEG for the same 
structure. Intermediate slices and the back surface are not shown. 
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Problem Structure 
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Figure 8. Front race 



856 




m 

o 

oo 



41 



Figure 



r~- 

co 




ki 



Figure 10. Freer Face o£ Structure 
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Back Face of Structure 



Table I 

Listing of Input Data Cards 
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Figure 12. Generated Mesh For Front Face 
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Figure 13. Generated Mesh For First Intermediate Slice 
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Figure 14. Generated Mesh For Second Intermediate Slice 
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Figure 15. Generated Mesh For Third Intermediate Slice 
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Figure 16. Generated Mesh For Back Face 
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Figure 17. Generated Mesh For Front Face 
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APPENDIX B 

A CLASSICAL PROBLEM SOLVED WITH TRISOP 



TRISOP is a relatively new computer system which uses isoparametric 
finite elements to obtain stresses, strains, and displacements of 
three-dimensional structures that are subjected to static loads. 

Although a number of test problems have been conducted with 
TRISOP, the input preparation was too voluminous to contemplate 
a full evaluation of the system. Without QUAMEG , the difficulties 
discussed in Chapter 1 limited the size and number of problems that 
could be solved. In order to evaluate TRISOP as a complete system, 
a classical problem was selected and a convergence study was made, 
using QUAMEG and TRISOP. 

The problem selected was that of the pinched, thin-shell cylinder. 
Obviously, this problem should be solved using thin-shell elements 
for maximum efficiency. The purpose of this test was to determine 
if the three-dimensional solution would converge to the same results 
as are reported in Ref. 3 and also to determine how well thin-shell 
action is captured by TRISOP. The results reported by Cantin [3] 
are included in Table II with the results obtained in' the test. The 
problem is described in Figure 18. 

Conclusions : 

Although the results of this convergence study indicate a 
tendency to approach the same results obtained by thin-shell theory, 
acceptable results can be obtained only with a very fine mesh. 
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TRISOP can be used to solve problems for which there are no 
closed-form or accurate numerical methods of solution available. 
However, the user must be prepared to pay the price of long compu- 
tational times, for some problems. 

Table II 

Displacement at the applied load for a pinched cylinder. 



Bogner et al. Cantin 



No. 

of 
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Mesh 


Displace. , 
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Figure 18. Pinched Cylinder 



APPENDIX C 

COMPUTER LISTING (QUAMEG) 
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APPENDIX D 

COMPUTER LISTING (CUMEG) 
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